

***********************************************************
* (1.0) Set up the survival analysis
***********************************************************			
	putexcel set  `"$d_tab\Table_A10.xlsx"', replace 
	putexcel C1 = "First Births", hcenter
	putexcel D1:E1 = "Second Births", merge hcenter
	putexcel F1:G1 = "Third and Higher Order Births", merge hcenter
	putexcel D2 = "Son",  hcenter
	putexcel E2 = "No Son",  hcenter
	putexcel F2 = "Son",  hcenter
	putexcel G2 = "No Son",  hcenter
	putexcel A2:G2, border(bottom)

	use 	"$f_2PKbirthpanel", clear	

	keep if rural == 1
	gen svydate		= (88*12)+7
				
	merge m:1 prov year using  "$d_data\IMR.dta", gen(IMRmerge)

	gen 	birth_year = .
	replace birth_year = year
		
	cap drop 	birth_month
	gen 		birth_month = .
	replace 	birth_month = pregmonth 

	replace 	marriage_year = 1900 + year_marriage

	gen 		marriage_age = age_marriage
	
	* Keeping only live births
	keep if pregresult<3
	
	* Parity is the order of live births, twins or multiples will be given a random parity
	* we'll use this method so that younger siblings of multiples are assigned appropriate parity
	* then we will drop the twin/multiples in the next step. 
	gen parity = order_livebirth - 1
	duplicates tag sample1pk id year pregmonth, gen(twins)
	replace mult = 1 if twins>0 & parity !=.	
	drop if mult == 1
	
	cap drop _merge
	merge 	m:1 prov year using `"$d_data\ProductionData.dta"'
	
	cap drop _merge
	merge m:1 prov using `"$d_data\sentdown_mean.dta"'
	drop _merge 
	
	gen 	culturalrev = 0 
	replace culturalrev = mean_sentdown if year >= 1966 & year <=1969

	keep  id parity cmale son svydate year_marriage month_marriage marcmc birth_year birth_month breastfeed_month LLFyear prov m_year_birth year educ marriage_age U5MR_5yr_ave grain agric_production gdp culturalrev pop_total_agric child_death_age

	gen eventyear = birth_year - LLFyear

	gen years_since_LLF = LLFyear - year 
	drop if years_since_LLF < -9 | years_since_LLF > 7
	drop if year >= 1980 | year < 1964

	replace breastfeed_month = . if breastfeed_month>96
	
	gen BF6 = 0 if breastfeed_month != . 
		replace BF6 = 1 if breastfeed_month >= 6 
	gen BF12 = 0 if breastfeed_month != . 
		replace BF12 = 1 if breastfeed_month >= 12 
	gen mort6 = child_death_age<6
	gen mort12 = child_death_age<12
		
	* Total births occurring and surviving to 6 months during LLF period 
	putexcel B3 = "Number of Births"
	* From vital statistics 
	local tot_births = 252763550 
	
	* We estimate the number of babies surviving to 6 months in each parity
	* and sex composition group by weighting for the share in each parity/sex comp group
	* and accounting for 0-5 month mortality rate 
	
		* First births
		count if cmale!=. 
			local num = r(N)
		count if parity==0
			local wt = r(N)/`num'
		sum mort6 if parity==0 
			local births_1 = round(`tot_births'*(1-r(mean))*`wt', 1)
			putexcel C3 = "`births_1'"
		
		* Second births, with older son 
		count if parity==1 & son==1
			local wt = r(N)/`num'
		sum mort6 if parity==1 & son==1
			local births_2s = round(`tot_births'*(1-r(mean))*`wt', 1)
			putexcel D3 = "`births_2s'"
			
		* Second births, with older son 
		count if parity==1 & son==0
			local wt = r(N)/`num'
		sum mort6 if parity==1 & son==0
			local births_2ns = round(`tot_births'*(1-r(mean))*`wt', 1)
			putexcel E3 = "`births_2ns'"

		* Third and higher births, with older son 
		count if parity>=2 & son==1
			local wt = r(N)/`num'
		sum mort6 if parity>=2 & son==1
			local births_3s = round(`tot_births'*(1-r(mean))*`wt', 1)
			putexcel F3 = "`births_3s'"
			
		* Third and higher births, with no older son 
		count if parity>=2 & son==0
			local wt = r(N)/`num'
		sum mort6 if parity>=2 & son==0
			local births_3ns = round(`tot_births'*(1-r(mean))*`wt', 1) 
			putexcel G3 = "`births_3ns'"
			
	* Overall share breastfed 
	putexcel A5 = "Breastfed Mortality Rate"
	putexcel B6 = "Overall Share Breastfed"
		sum BF12 if parity==0  & years_since_LLF<0 & mort12==0
			local BF12_1 = round(r(mean), .001)
			scalar cell = string(`BF12_1', "%05.3f") 
			putexcel C6 = "`=cell'", hcenter
		sum BF12 if parity==1 & son==1 & years_since_LLF<0 & mort12==0
			local BF12_2s = round(r(mean), .001)
			scalar cell = string(`BF12_2s', "%05.3f") 
			putexcel D6 = "`=cell'", hcenter
		sum BF12 if parity==1 & son==0 & years_since_LLF<0 & mort12==0
			local BF12_2ns = round(r(mean), .001)
			scalar cell = string(`BF12_2ns', "%05.3f") 
			putexcel E6 = "`=cell'", hcenter
		sum BF12 if parity>=2 & son==1 & years_since_LLF<0 & mort12==0
			local BF12_3s = round(r(mean), .001)
			scalar cell = string(`BF12_3s', "%05.3f") 
			putexcel F6 = "`=cell'", hcenter
		sum BF12 if parity>=2 & son==0 & years_since_LLF<0 & mort12==0
			local BF12_3ns = round(r(mean), .001)
			scalar cell = string(`BF12_3ns', "%05.3f") 
			putexcel G6 = "`=cell'", hcenter
			
	* Mortality rate between 6-11 
	putexcel B7 = "6-12 Month Mortality Rate "
	sum mort12 if mort6==0 & parity==0 &  years_since_LLF<0
			local mort_1 = round(r(mean), .001)
			scalar cell = string(`mort_1', "%05.3f") 
			putexcel C7 = "`=cell'", hcenter
	sum mort12 if mort6==0 & parity==1 & son==1 & years_since_LLF<0
			local mort_2s = round(r(mean), .001)
			scalar cell = string(`mort_2s', "%05.3f") 
			putexcel D7 = "`=cell'", hcenter
	sum mort12 if mort6==0 & parity==1 & son==0 & years_since_LLF<0
			local mort_2ns = round(r(mean), .001)
			scalar cell = string(`mort_2ns', "%05.3f") 
			putexcel E7 = "`=cell'", hcenter
	sum mort12 if mort6==0 & parity>=2 & son==1 & years_since_LLF<0
			local mort_3s = round(r(mean), .001)
			scalar cell = string(`mort_3s', "%05.3f") 
			putexcel F7 = "`=cell'", hcenter
	sum mort12 if mort6==0 & parity>=2 & son==0 & years_since_LLF<0
			local mort_3ns = round(r(mean), .001)
			scalar cell = string(`mort_3ns', "%05.3f") 
			putexcel G7 = "`=cell'", hcenter
	
	* Breastfed mortality rate between 6-11 
	putexcel B9 = "Breastfed 6-12 Month Mortality Rate "
	local BFmort_1  = round(`mort_1'/(`BF12_1' + 1.8*(1-`BF12_1')), .001)
			scalar cell = string(`BFmort_1', "%05.3f") 
			putexcel C9 = "`=cell'", hcenter
	local BFmort_2s  = round(`mort_2s'/(`BF12_2s' + 1.8*(1-`BF12_2s')), .001)
			scalar cell = string(`BFmort_2s', "%05.3f") 
			putexcel D9 = "`=cell'", hcenter
	local BFmort_2ns  = round(`mort_2ns'/(`BF12_2ns' + 1.8*(1-`BF12_2ns')), .001)
			scalar cell = string(`BFmort_2s', "%05.3f") 
			putexcel E9 = "`=cell'", hcenter
	local BFmort_3s  = round(`mort_3s'/(`BF12_3s' + 1.8*(1-`BF12_3s')), .001)
			scalar cell = string(`BFmort_3s', "%05.3f") 
			putexcel D9 = "`=cell'", hcenter
	local BFmort_3ns  = round(`mort_3ns'/(`BF12_3ns' + 1.8*(1-`BF12_3ns')), .001)
			scalar cell = string(`BFmort_3s', "%05.3f") 
			putexcel E9 = "`=cell'", hcenter
			
	* Non breastfed mortality rate between 6-11 
	putexcel B10 = "Non-Breastfed 6-12 Month Mortality Rate "
	local nBFmort_1  = round(`BFmort_1'*1.8, .001)
			scalar cell = string(`nBFmort_1', "%05.3f") 
			putexcel C10 = "`=cell'", hcenter
	local nBFmort_2s  = round(`BFmort_2s'*1.8, .001)
			scalar cell = string(`nBFmort_2s', "%05.3f") 
			putexcel D10 = "`=cell'", hcenter
	local nBFmort_2ns  = round(`BFmort_2ns'*1.8, .001)
			scalar cell = string(`nBFmort_2ns', "%05.3f") 
			putexcel E10 = "`=cell'", hcenter
	local nBFmort_3s  =  round(`BFmort_3s'*1.8, .001)
			scalar cell = string(`nBFmort_3s', "%05.3f") 
			putexcel F10 = "`=cell'", hcenter
	local nBFmort_3ns  = round(`BFmort_3ns'*1.8, .001)
			scalar cell = string(`nBFmort_3s', "%05.3f") 
			putexcel G10 = "`=cell'", hcenter
	
	* Breastfeeding Gap 
	putexcel B10:G10, border(bottom)
	putexcel B11 = "Breastfeeding Mortality Delta"
	local mortGap_1  = round(`nBFmort_1' -`BFmort_1', .001)
			scalar cell = string(`mortGap_1', "%05.3f") 
			putexcel C11 = "`=cell'", hcenter
	local mortGap_2s  = round(`nBFmort_2s' -`BFmort_2s', .001)
			scalar cell = string(`mortGap_2s', "%05.3f") 
			putexcel D11 = "`=cell'", hcenter
	local mortGap_2ns  = round(`nBFmort_2ns' -`BFmort_2ns', .001)
			scalar cell = string(`mortGap_2ns', "%05.3f") 
			putexcel E11 = "`=cell'", hcenter
	local mortGap_3s  = round(`nBFmort_3ns' -`BFmort_3s', .001)
			scalar cell = string(`mortGap_3s', "%05.3f") 
			putexcel F11 = "`=cell'", hcenter			
	local mortGap_3ns  = round(`nBFmort_3ns' -`BFmort_3ns', .001)
			scalar cell = string(`mortGap_3ns', "%05.3f") 
			putexcel G11 = "`=cell'", hcenter
	
	* Gender gap in breastfeeding
	putexcel A13 = "Gender Gap in Breastfeeding"
	putexcel B14 = "Gender Gap in Breastfeeding, Pre LLF Period"
	qui sum BF12 if years_since_LLF<0 & parity==0 & cmale==1 
		local male = r(mean)
	qui sum BF12 if years_since_LLF<0 & parity==0 & cmale==0 
		local female = r(mean)
		local BFgap_1 = round(`male'-`female', .001)
		scalar cell = string(`BFgap_1', "%05.3f")
		putexcel C14 = "`=cell'", hcenter
		
	qui sum BF12 if years_since_LLF<0 & parity==1 & cmale==1 & son==1
		local male = r(mean)
	qui sum BF12 if years_since_LLF<0 & parity==1 & cmale==0 & son==1
		local female = r(mean)
		local BFgap_2s = round(`male'-`female', .001)
		scalar cell = string(`BFgap_2s', "%05.3f")
		putexcel D14 = "`=cell'", hcenter
		
	qui sum BF12 if years_since_LLF<0 & parity==1 & cmale==1 & son==0
		local male = r(mean)
	qui sum BF12 if years_since_LLF<0 & parity==1 & cmale==0 & son==0
		local female = r(mean)
		local BFgap_2ns = round(`male'-`female', .001)
		scalar cell = string(`BFgap_2ns', "%05.3f")
		putexcel E14 = "`=cell'", hcenter

	qui sum BF12 if years_since_LLF<0 & parity>=2 & cmale==1 & son==1
		local male = r(mean)
	qui sum BF12 if years_since_LLF<0 & parity>=2 & cmale==0 & son==1
		local female = r(mean)
		local BFgap_3s = round(`male'-`female', .001)
		scalar cell = string(`BFgap_3s', "%05.3f")
		putexcel F14 = "`=cell'", hcenter
		
	qui sum BF12 if years_since_LLF<0 & parity>=2 & cmale==1 & son==0
		local male = r(mean)
	qui sum BF12 if years_since_LLF<0 & parity>=2 & cmale==0 & son==0
		local female = r(mean)
		local BFgap_3ns = round(`male'-`female', .001)
		scalar cell = string(`BFgap_3ns', "%05.3f")
		putexcel G14 = "`=cell'", hcenter
		
	* Parameters from sex composition strategy use analysis 
	putexcel B15 = "Sex Composition Strategy Use, Pre LLF Period"
		putexcel C15 = "0.003", hcenter
		putexcel D15 = "0.012", hcenter
		putexcel E15 = "0.054", hcenter
		putexcel F15 = "0.053", hcenter
		putexcel G15 = "0.110", hcenter

	putexcel B16 = "Gender Gap in Breastfeeding / Sex Composition Strategy Use"
		local preratio_1 = `BFgap_1'/0.003
		scalar cell = string(`preratio_1', "%05.3f")
		putexcel C16 = "`=cell'", hcenter
		
		local preratio_2s = `BFgap_2s'/0.012
		scalar cell = string(`preratio_2s', "%05.3f")
		putexcel D16 = "`=cell'", hcenter
		
		local preratio_2ns = `BFgap_2ns'/0.054
		scalar cell = string(`preratio_2ns', "%05.3f")
		putexcel E16 = "`=cell'", hcenter
		
		local preratio_3s = `BFgap_3s'/0.012
		scalar cell = string(`preratio_3s', "%05.3f")
		putexcel F16 = "`=cell'", hcenter
		
		local preratio_3ns = `BFgap_3ns'/0.054
		scalar cell = string(`preratio_3ns', "%05.3f")
		putexcel G16 = "`=cell'", hcenter

	* Parameters from sex composition strategy use analysis 
	putexcel B18 = "Sex Composition Strategy Use, Post LLF Period"
		putexcel C18 = "0.100", hcenter
		putexcel D18 = "0.067", hcenter
		putexcel E18 = "0.271", hcenter
		putexcel F18 = "0.065", hcenter
		putexcel G18 = "0.330", hcenter
	putexcel B18:G18, border(bottom)
	
	putexcel B20 = "Reported Gender Gap in Breastfeeding, Late LLF Period"
		qui sum BF12 if years_since_LLF>=0 & parity==0 & cmale==1 
			local male = r(mean)
		qui sum BF12 if years_since_LLF>=0 & parity==0 & cmale==0 
			local female = r(mean)
		local lBFgap_1 = round(`male'-`female', .001)
		scalar cell = string(`lBFgap_1', "%05.3f")
		putexcel C20 = "`=cell'", hcenter

		qui sum BF12 if years_since_LLF>=0 & parity==1 & cmale==1 & son==1
			local male = r(mean)
		qui sum BF12 if years_since_LLF>=0 & parity==1 & cmale==0 & son==1
			local female = r(mean)
		local lBFgap_2s = round(`male'-`female', .001)
		scalar cell = string(`lBFgap_2s', "%05.3f")
		putexcel D20 = "`=cell'", hcenter

		qui sum BF12 if years_since_LLF>=0 & parity==1 & cmale==1 & son==0
			local male = r(mean)
		qui sum BF12 if years_since_LLF>=0 & parity==1 & cmale==0 & son==0
			local female = r(mean)
		local lBFgap_2ns = round(`male'-`female', .001)
		scalar cell = string(`lBFgap_2ns', "%05.3f")
		putexcel E20 = "`=cell'", hcenter

		qui sum BF12 if years_since_LLF>=0 & parity>=2 & cmale==1 & son==1
			local male = r(mean)
		qui sum BF12 if years_since_LLF>=0 & parity>=2 & cmale==0 & son==1
			local female = r(mean)
		local lBFgap_3s = round(`male'-`female', .001)
		scalar cell = string(`lBFgap_3s', "%05.3f")
		putexcel F20 = "`=cell'", hcenter

		qui sum BF12 if years_since_LLF>=0 & parity>=2 & cmale==1 & son==0
			local male = r(mean)
		qui sum BF12 if years_since_LLF>=0 & parity>=2 & cmale==0 & son==0
			local female = r(mean)
		local lBFgap_3ns = round(`male'-`female', .001)
		scalar cell = string(`lBFgap_3ns', "%05.3f")
		putexcel G20 = "`=cell'", hcenter
		
	putexcel A22 = "Missing Girls Due to Breastfeeding"
	putexcel B26 = "Breastfeeding Mortality Gap X Reported Gender Gap in Breastfeeding, Post LLF Period"
		local cell_1 = `mortGap_1'*`lBFgap_1'
		scalar cell_1 = string(`cell_1', "%05.3f")
		putexcel C26 = "`=cell_1'", hcenter
		
		local cell_2s = `mortGap_2s'*`lBFgap_2s'
		scalar cell = string(`cell_2s', "%05.3f")
		putexcel D26 = "`=cell'", hcenter
		
		local cell_2ns = `mortGap_2ns'*`lBFgap_2ns'
		scalar cell = string(`cell_2ns', "%05.3f")
		putexcel E26 = "`=cell'", hcenter
		
		local cell_3s = `mortGap_3s'*`lBFgap_3s'
		scalar cell = string(`cell_3s', "%05.3f")
		putexcel F26 = "`=cell'", hcenter
		
		local cell_3ns = `mortGap_3ns'*`lBFgap_3ns'
		scalar cell = string(`cell_3ns', "%05.3f")
		putexcel G26 = "`=cell'", hcenter
		
	putexcel B27 = "Excess Female Deaths Due to Reported Breastfeeing Gap"
		local excess_1 = round(`cell_1'*`births_1', 1)
		putexcel C27 = "`excess_1'", hcenter
		local excess_2s = round(`cell_2s'*`births_2s', 1)
		putexcel D27 = "`excess_2s'", hcenter
		local excess_2ns = round(`cell_2ns'*`births_2ns', 1)
		putexcel E27 = "`excess_2ns'", hcenter
		local excess_3s = round(`cell_3s'*`births_3s', 1)
		putexcel F27 = "`excess_3s'", hcenter
		local excess_3ns = round(`cell_3ns'*`births_3ns', 1)
		putexcel G27 = "`excess_3ns'", hcenter	
		
	putexcel A27:G27, border(bottom)
	
**********************************************************
* (2.0) Discrete time hazard model of breastfeeding duration 
**********************************************************

{	
	do `"$d_do\LLF_StackData.do"'

	gen years_since_LLF = LLFyear - year 
	drop if years_since_LLF < -9 | years_since_LLF > 7
	drop if year >=1980
	drop if year <1964

	cap drop uniqueid 
	egen uniqueid = group(StackDataID id)
	replace year_marriage = 1900 + year_marriage

	cap drop years_since_LLF
	gen years_since_LLF = year - LLFyear

	global controls1 	"i.educ U5MR_5yr_ave"
	global controls2 	"grain agric_production gdp pop_total_agric culturalrev"
	
	*********************
	* (2.1) First births
	*********************
	preserve
	
	keep if parity == 0  & year_marriage >= 1960 & birth_year <= 1979
	keep marriage_age cmale m_year_birth educ U5MR_5yr_ave $controls2 birth_month birth_year breastfeed_month LLFyear year uniqueid prov svydate marcmc FE_t FE_prov
	
	* Make this a woman-month level dataset, allowing LLF to vary over the interval
	* Gen monthly indicators to capture survival time - equal to zero until the breastfeeding interval ends, 1 thereafter
	gen 	month0 = 1 if breastfeed_month == 1 & breastfeed_month != .
	
	forval i = 1/60 { 
		* The indictor will be zero for all months before breastfeeding stops and 1 the month breastfeeding ends
		gen 	month`i'=0 if breastfeed_month > `i' & breastfeed_month != .
		replace month`i'=1 if breastfeed_month == `i'  	
	}
	
	replace month60 = 1 if breastfeed_month > 60 & breastfeed_month != .
	
	* Identify cmc time for each month in the interval to determine if LLF was active
	gen start_cmc = 12*(birth_year-1900) + birth_month
	
	forval i = 1/60 { 
		cap drop cmc`i'
		gen cmc`i' = .
		
		replace cmc`i' = start_cmc + `i' if month`i' != .
		
		cap drop pd_year`i'  
		gen pd_year`i' = . 
		replace pd_year`i' = cmc`i'/12 + 1900 if month`i' != .
		
		cap drop pd_LLF`i' 
		gen pd_LLF`i' = . 
		replace pd_LLF`i' = 0 if month`i' != . 
		replace pd_LLF`i' = 1 if pd_year`i' >= LLFyear & month`i' != . 
		} 
		
	drop cmc* start_cmc
	
	reshape long pd_LLF pd_year month, i(uniqueid) j(interval_month)
	replace pd_year = floor(pd_year)
	rename month end_bf
	rename interval_month month
	stset month, failure(end_bf) id(uniqueid)
	drop if end_bf==.
	
	tab month, gen(month_)
			
		* Generate 2 way period * LLF interactions 
		forval i = 1/60 { 
			gen LLF_month_`i' = pd_LLF*month_`i' 
			} 

		* Generate 2 way period * male interaction
		forval i = 1/60 { 
			gen month_`i'_male = month_`i'*cmale 
			} 		
		
		* Generate 2 way LLF * male interaction
		gen pd_LLF_male = pd_LLF*cmale
		
		* Generate 3 way period * LLF * male interactions 
		forval i = 1/60 { 
			gen LLF_month_`i'_male = pd_LLF*month_`i'*cmale 
			} 
			 	
	gen birth_year2LLF = birth_year-LLFyear
	
	keep if pd_year >= 1964
	keep if birth_year2LLF >= -4	
	keep if birth_year2LLF <= 9		
	
	xi: logit 	end_bf 	pd_LLF cmale month_1-month_11 month_13-month_60 ///
						LLF_month_1-LLF_month_11 LLF_month_13-LLF_month_60 pd_LLF_male ///
						month_1_male-month_11_male month_13_male-month_60_male /// 
						LLF_month_1_male-LLF_month_11_male LLF_month_13_male-LLF_month_60_male ///
						m_year_birth $controls1 $controls2 ///
						i.prov i.pd_year if pd_year!=1967 & pd_year != 1968, cluster(prov) ltolerance(1e-1) iter(50)	

	* Populate a Life Table using point estimates, using preLLF as the reference group
	keep if e(sample)
	keep if pd_LLF==0
	
		* The reference group is males before LLF, 
		*  set all LLF indicators and interactions to 0, 
		*  all male indicators to 1
	
		* Replace LLF = 0 
		replace pd_LLF = 0

		* Replace male indicators to 1 
		replace cmale = 1 

		* Replace LLF*month interactions with zero and predict reference group
		forval i = 1/60 { 			
			replace LLF_month_`i' = 0
			} 
			
		* Replace LLF*male interacations with zero 
		replace pd_LLF_male = 0 
		
		* Generate new month*male indicators using cmale == 1 
		forval i = 1/60 { 			
			replace month_`i'_male = month_`i'*cmale
			}

		* Replace LLF*male*month interacations with zero 
		forval i = 1/60 { 			
			replace LLF_month_`i'_male = 0
			}
		
		* Predict the outcome in the reference group
		cap drop xb
		predict xb, xb
		
		cap drop qx
		predict qx
		
			forval q = 1/60 { 
			
				sum xb if month == `q' 
				local xb = r(mean)
				
				sum qx if month == `q' 
				local qx = r(mean)
				
				if `q' == 1 {
					mat LT_NoLLF_male = `q', `xb', `qx' 
					} 
				else if `q' > 1 { 
					mat temp = `q', `xb', `qx' 
					mat LT_NoLLF_male = LT_NoLLF_male \ temp
					}
				}
				
		* FEMALES PRE LLF
		* Keep all LLF indicators and interactions to 0, 
		* all male indicators now to 0
	
		* Replace male indicators to 0 
		replace cmale = 0 
		
		* Generate new month*male indicators using cmale == 1 
		forval i = 1/60 { 			
			replace month_`i'_male = 0
			}
		
		* Now predict the outcome
		cap drop xb
		predict xb, xb
		
		cap drop qx
		predict qx
		
			forval q = 1/60 { 
			
				sum xb if month == `q' 
				local xb = r(mean)
				
				sum qx if month == `q' 
				local qx = r(mean)
				
				if `q' == 1 {
					mat LT_NoLLF_female = `q', `xb', `qx' 
					} 
				else if `q' > 1 { 
					mat temp = `q', `xb', `qx' 
					mat LT_NoLLF_female = LT_NoLLF_female \ temp
					}
				}

		* MALES POST LLF
		* Set all LLF indicators and interactions to 1, 
		* all male indicators to 1
	
		* Replace LLF = 1 
		replace pd_LLF = 1

		* Replace male indicator with 1
		replace cmale = 1
		
		* Remake LLF*month interactions with pd_LLF==1
		forval i = 1/60 { 			
			replace LLF_month_`i' = month_`i'*pd_LLF 
			} 
			
		* Remake LLF*male interacations with 1 
		replace pd_LLF_male = 1 
		
		* Generate new month*male indicators using cmale == 1 
		forval i = 1/60 { 			
			replace month_`i'_male = month_`i'*cmale
			}
			
		* Remake LLF*male*month interacations with cmale == 1 and LLF == 1
		forval i = 1/60 { 			
			replace LLF_month_`i'_male = month_`i'*cmale*pd_LLF
			}
			
		* Predict the outcome 
		cap drop xb
		predict xb, xb
		
		cap drop qx
		predict qx
		
			forval q = 1/60 { 
			
				sum xb if month == `q' 
				local xb = r(mean)
				
				sum qx if month == `q' 
				local qx = r(mean)
				
				if `q' == 1 {
					mat LT_LLF_male = `q', `xb', `qx' 
					} 
				else if `q' > 1 { 
					mat temp = `q', `xb', `qx' 
					mat LT_LLF_male = LT_LLF_male \ temp
					}
				}

		
		* FEMALES POST LLF
		* Set all LLF indicators and interactions with month to 1, 
		* all male indicators to 0

		* Replace male indicator with 0
		replace cmale = 0
		
		* Replace LLF*male interacations with 0 
		replace pd_LLF_male = 0
		
		* Replace month*male indicators with 0 
		forval i = 1/60 { 			
			replace month_`i'_male = 0
			}
			
		* Replace LLF*male*month interacations to 0
		forval i = 1/60 { 			
			replace LLF_month_`i'_male = 0
			}
			
		* Predict the outcome 
		cap drop xb
		predict xb, xb
		
		cap drop qx
		predict qx
		
			forval q = 1/60 { 
			
				sum xb if month == `q' 
				local xb = r(mean)
				
				sum qx if month == `q' 
				local qx = r(mean)
				
				if `q' == 1 {
					mat LT_LLF_female = `q', `xb', `qx' 
					} 
				else if `q' > 1 { 
					mat temp = `q', `xb', `qx' 
					mat LT_LLF_female = LT_LLF_female \ temp
					}
				}				

	* Compute Lx and Dx for each month, males and females, with and without LLF
	mat temp = J(60, 1, .) 
	mat LT_NoLLF_male 	= LT_NoLLF_male, temp, temp
	mat LT_NoLLF_female = LT_NoLLF_female, temp, temp
	mat LT_LLF_male 	= LT_LLF_male, temp, temp
	mat LT_LLF_female 	= LT_LLF_female, temp, temp
		
	mat LT_NoLLF_male[1,4] 		= 1
	mat LT_NoLLF_female[1,4] 	= 1
	mat LT_LLF_male[1,4] 		= 1
	mat LT_LLF_female[1,4] 		= 1
	
	forval i = 1/60 {
		mat LT_NoLLF_male[`i',5] 	= LT_NoLLF_male[`i',3]*LT_NoLLF_male[`i',4]
		mat LT_NoLLF_female[`i',5] 	= LT_NoLLF_female[`i',3]*LT_NoLLF_female[`i',4]
		mat LT_LLF_male[`i',5] 		= LT_LLF_male[`i',3]*LT_LLF_male[`i',4]
		mat LT_LLF_female[`i',5] 	= LT_LLF_female[`i',3]*LT_LLF_female[`i',4]
		
		if `i' < 60 { 
			local j = `i' + 1
			mat LT_NoLLF_male[`j' , 4] 		= LT_NoLLF_male[`i' , 4] - LT_NoLLF_male[`i',5]
			mat LT_NoLLF_female[`j' , 4] 	= LT_NoLLF_female[`i' , 4] - LT_NoLLF_female[`i',5]
			mat LT_LLF_male[`j' , 4] 		= LT_LLF_male[`i' , 4] - LT_LLF_male[`i',5]
			mat LT_LLF_female[`j' , 4] 		= LT_LLF_female[`i' , 4] - LT_LLF_female[`i',5]
		} 
	}
		
	mat colnames LT_NoLLF_male 		= "Month" "XB_NoLLF_male" "Qx_NoLLF_male" "Lx_NoLLF_male" "Dx_NoLLF_male" 
	mat colnames LT_NoLLF_female 	= "Month" "XB_NoLLF_female" "Qx_NoLLF_female" "Lx_NoLLF_female" "Dx_NoLLF_female" 
	mat colnames LT_LLF_male 		= "Month" "XB_LLF_male" "Qx_LLF_male" "Lx_LLF_male" "Dx_LLF_male"
	mat colnames LT_LLF_female 		= "Month" "XB_LLF_female" "Qx_LLF_female" "Lx_LLF_female" "Dx_LLF_female"
	
	mat temp1 = LT_NoLLF_female[1..., 2...]	
	mat temp2 = LT_LLF_male[1..., 2...]	
	mat temp3 = LT_LLF_female[1..., 2...]
	
	mat colnames temp1 	= "XB_NoLLF_female" "Qx_NoLLF_female" "Lx_NoLLF_female" "Dx_NoLLF_female" 
	mat colnames temp2 	= "XB_LLF_male" "Qx_LLF_male" "Lx_LLF_male" "Dx_LLF_male"
	mat colnames temp3 	= "XB_LLF_female" "Qx_LLF_female" "Lx_LLF_female" "Dx_LLF_female"
	
	mat LifeTable_1 = LT_NoLLF_male, temp1, temp2, temp3
	
	* Plot the CDF
	clear 
	svmat LifeTable_1, names(col)
	
	local mean_NoLLF_male = 0
	local dx_cum_NoLLF_male = 0 
	local flag_NoLLF_male = 1

	forval i = 1/26 { 
		local j = `i' - 1
		local int = LT_NoLLF_male[`i', 1]*LT_NoLLF_male[`i', 5]
		local dx = LT_NoLLF_male[`i', 5]
		local mean_NoLLF_male = `mean_NoLLF_male' + `int' 
		local dx_cum_NoLLF_male = `dx_cum_NoLLF_male' + `dx' 
		
		local y = LT_NoLLF_male[`i', 4]
		if `y' < .5 & `flag_NoLLF_male' == 1 { 
			local flag_NoLLF_male = 2 
			local y_1 = LT_NoLLF_male[`j', 4]
			local x = LT_NoLLF_male[`i', 1]
			local beta = `y' - `y_1'
			local alpha = `y'-(`beta'*`x')
			local median_NoLLF_male = (.5-`alpha')/`beta'
			}
	}		

	local mean_NoLLF_female = 0
	local dx_cum_NoLLF_female = 0 
	local flag_NoLLF_female = 1

	forval i = 1/26 { 
		local j = `i' - 1
		local int = LT_NoLLF_female[`i', 1]*LT_NoLLF_female[`i', 5]
		local dx = LT_NoLLF_female[`i', 5]
		local mean_NoLLF_female = `mean_NoLLF_female' + `int' 
		local dx_cum_NoLLF_female = `dx_cum_NoLLF_female' + `dx' 
		
		local y = LT_NoLLF_female[`i', 4]
		if `y' < .5 & `flag_NoLLF_female' == 1 { 
			local flag_NoLLF_female = 2 
			local y_1 = LT_NoLLF_female[`j', 4]
			local x = LT_NoLLF_female[`i', 1]
			local beta = `y' - `y_1'
			local alpha = `y'-(`beta'*`x')
			local median_NoLLF_female = (.5-`alpha')/`beta'
			}
	}		
	
	local mean_LLF_male = 0 
	local dx_cum_LLF_male = 0 
	local flag_LLF_male = 1
	forval i = 1/28 { 
		local j = `i' - 1
		local int = LT_LLF_male[`i', 1]*LT_LLF_male[`i', 5]
		local dx = LT_LLF_male[`i', 5]
		local mean_LLF_male = `mean_LLF_male' + `int' 
		local dx_cum_LLF_male = `dx_cum_LLF_male' + `dx' 
		
		local y = LT_LLF_male[`i', 4]
		if `y' < .5 & `flag_LLF_male' == 1 { 
			local flag_LLF_male = 2 
			local y_1 = LT_LLF_male[`j', 4]
			local x = LT_LLF_male[`i', 1]
			local beta = `y' - `y_1'
			local alpha = `y'-(`beta'*`x')
			local median_LLF_male = (.5-`alpha')/`beta'
			}
		} 

	local mean_LLF_female = 0 
	local dx_cum_LLF_female = 0 
	local flag_LLF_female = 1
	forval i = 1/28 { 
		local j = `i' - 1
		local int = LT_LLF_female[`i', 1]*LT_LLF_female[`i', 5]
		local dx = LT_LLF_female[`i', 5]
		local mean_LLF_female = `mean_LLF_female' + `int' 
		local dx_cum_LLF_female = `dx_cum_LLF_female' + `dx' 
		
		local y = LT_LLF_female[`i', 4]
		if `y' < .5 & `flag_LLF_female' == 1 { 
			local flag_LLF_female = 2 
			local y_1 = LT_LLF_female[`j', 4]
			local x = LT_LLF_female[`i', 1]
			local beta = `y' - `y_1'
			local alpha = `y'-(`beta'*`x')
			local median_LLF_female = (.5-`alpha')/`beta'
			}
		} 
		
	local mean_NoLLF_male 	= `mean_NoLLF_male'*`dx_cum_NoLLF_male'
	local mean_NoLLF_female = `mean_NoLLF_female'*`dx_cum_NoLLF_female'
	local mean_LLF_male 	= `mean_LLF_male'*`dx_cum_LLF_male'
	local mean_LLF_female 	= `mean_LLF_female'*`dx_cum_LLF_female'
	
	local change_in_gap = (`mean_LLF_male'-`mean_LLF_female')-(`mean_NoLLF_male'-`mean_NoLLF_female')
	di "`mean_NoLLF_male' | `mean_NoLLF_female' | `mean_LLF_male' | `mean_LLF_female'"
	mat BF_effects_First = `change_in_gap'
	
	restore
	
	****************************************
	* (2.2) Second births, No Previous Sons
	****************************************
			
	preserve
		
	keep if parity == 1 & year_marriage >= 1960 & birth_year <= 1979 & son == 0
	keep marriage_age cmale m_year_birth educ U5MR_5yr_ave $controls2 birth_month birth_year breastfeed_month LLFyear year uniqueid prov svydate marcmc
		
	gen 	month0 = 1 if breastfeed_month == 1 & breastfeed_month != .
		
		forval i = 1/60 { 
			gen 	month`i' = 0 if breastfeed_month > `i' & breastfeed_month != .
			replace month`i' = 1 if breastfeed_month == `i'  	
		}
		
	replace month60 = 1 if breastfeed_month > 60 & breastfeed_month != .
		
	gen start_cmc = 12*(birth_year-1900) + birth_month
		
		forval i = 1/60 { 
			cap drop cmc`i'
			gen cmc`i' = .
			
			replace cmc`i' = start_cmc + `i' if month`i' != .
			
			cap drop pd_year`i'  
			gen pd_year`i' = . 
			replace pd_year`i' = cmc`i'/12 + 1900 if month`i' != .
			
			cap drop pd_LLF`i' 
			gen pd_LLF`i' = . 
			replace pd_LLF`i' = 0 if month`i' != . 
			replace pd_LLF`i' = 1 if pd_year`i' >= LLFyear & month`i' != . 
			} 
			
	drop cmc* start_cmc
		
	reshape long pd_LLF pd_year month, i(uniqueid) j(interval_month)
	replace pd_year = floor(pd_year)
	rename month end_bf
	rename interval_month month
	stset month, failure(end_bf) id(uniqueid)
	drop if end_bf==.
			
	tab month, gen(month_)
	forval i = 1/60 { 
		gen LLF_month_`i' = pd_LLF*month_`i' 
		} 

	forval i = 1/60 { 
		gen month_`i'_male = month_`i'*cmale 
		} 		
		
	gen pd_LLF_male = pd_LLF*cmale
		
	forval i = 1/60 { 
		gen LLF_month_`i'_male = pd_LLF*month_`i'*cmale 
		} 
			
	gen birth_year2LLF = birth_year-LLFyear
		
	keep if pd_year >= 1964
	keep if birth_year2LLF >= -4	
	keep if birth_year2LLF <= 9		
		
	xi: logit 	end_bf 	pd_LLF cmale month_1-month_11 month_13-month_60 ///
						LLF_month_1-LLF_month_11 LLF_month_13-LLF_month_60 pd_LLF_male ///
						month_1_male-month_11_male month_13_male-month_60_male /// 
						LLF_month_1_male-LLF_month_11_male LLF_month_13_male-LLF_month_60_male ///
						m_year_birth $controls1 $controls2 ///
						i.prov i.pd_year if pd_year!=1967 & pd_year != 1968, cluster(prov) ltolerance(1e-1) iter(30)	

	* REFERENCE GROUP: MALES PRE LLF
	replace pd_LLF = 0
	replace cmale = 1 

	forval i = 1/60 { 			
		replace LLF_month_`i' = 0
		} 
			
	replace pd_LLF_male = 0 
		
	forval i = 1/60 { 			
		replace month_`i'_male = month_`i'*cmale
		}

	forval i = 1/60 { 			
		replace LLF_month_`i'_male = 0
		}
		
	cap drop xb
	predict xb, xb
		
	cap drop qx
	predict qx
		
		forval q = 1/60 { 
			
			sum xb if month == `q' 
			local xb = r(mean)
				
			sum qx if month == `q' 
			local qx = r(mean)
				
			if `q' == 1 {
				mat LT_NoLLF_male = `q', `xb', `qx' 
				} 
			else if `q' > 1 { 
				mat temp = `q', `xb', `qx' 
				mat LT_NoLLF_male = LT_NoLLF_male \ temp
				}
			}
				
	* FEMALES PRE LLF
	replace cmale = 0 
		
	forval i = 1/60 { 			
		replace month_`i'_male = 0
		}
		
	cap drop xb
	predict xb, xb
		
	cap drop qx
	predict qx
		
		forval q = 1/60 { 
			
			sum xb if month == `q' 
			local xb = r(mean)
				
			sum qx if month == `q' 
			local qx = r(mean)
				
			if `q' == 1 {
				mat LT_NoLLF_female = `q', `xb', `qx' 
				} 
			else if `q' > 1 { 
				mat temp = `q', `xb', `qx' 
				mat LT_NoLLF_female = LT_NoLLF_female \ temp
				}
			}

	* MALES POST LLF
	replace pd_LLF = 1
	replace cmale = 1
		
	forval i = 1/60 { 			
		replace LLF_month_`i' = month_`i'*pd_LLF 
		} 
			
	replace pd_LLF_male = 1 
		
	forval i = 1/60 { 			
		replace month_`i'_male = month_`i'*cmale
		}
			
	forval i = 1/60 { 			
		replace LLF_month_`i'_male = month_`i'*cmale*pd_LLF
		}
			
	cap drop xb
	predict xb, xb
		
	cap drop qx
	predict qx
		
		forval q = 1/60 { 
			
			sum xb if month == `q' 
			local xb = r(mean)
				
			sum qx if month == `q' 
			local qx = r(mean)
				
			if `q' == 1 {
				mat LT_LLF_male = `q', `xb', `qx' 
				} 
			else if `q' > 1 { 
				mat temp = `q', `xb', `qx' 
				mat LT_LLF_male = LT_LLF_male \ temp
				}
			}

	* FEMALES POST LLF
	replace cmale = 0		
	replace pd_LLF_male = 0
		
	forval i = 1/60 { 			
		replace month_`i'_male = 0
		}
			
	forval i = 1/60 { 			
		replace LLF_month_`i'_male = 0
		}
			
	cap drop xb
	predict xb, xb
		
	cap drop qx
	predict qx
		
		forval q = 1/60 { 
			
			sum xb if month == `q' 
			local xb = r(mean)
				
			sum qx if month == `q' 
			local qx = r(mean)
				
			if `q' == 1 {
				mat LT_LLF_female = `q', `xb', `qx' 
				} 
			else if `q' > 1 { 
				mat temp = `q', `xb', `qx' 
				mat LT_LLF_female = LT_LLF_female \ temp
				}
			}			

	* Compute Lx and Dx for each month, males and females, with and without LLF
	mat temp = J(60, 1, .) 
	mat LT_NoLLF_male 	= LT_NoLLF_male, temp, temp
	mat LT_NoLLF_female = LT_NoLLF_female, temp, temp
	mat LT_LLF_male 	= LT_LLF_male, temp, temp
	mat LT_LLF_female 	= LT_LLF_female, temp, temp
		
	mat LT_NoLLF_male[1,4] 		= 1
	mat LT_NoLLF_female[1,4] 	= 1
	mat LT_LLF_male[1,4] 		= 1
	mat LT_LLF_female[1,4] 		= 1
	
	forval i = 1/60 {
		mat LT_NoLLF_male[`i',5] 	= LT_NoLLF_male[`i',3]*LT_NoLLF_male[`i',4]
		mat LT_NoLLF_female[`i',5] 	= LT_NoLLF_female[`i',3]*LT_NoLLF_female[`i',4]
		mat LT_LLF_male[`i',5] 		= LT_LLF_male[`i',3]*LT_LLF_male[`i',4]
		mat LT_LLF_female[`i',5] 	= LT_LLF_female[`i',3]*LT_LLF_female[`i',4]
		
		if `i' < 60 { 
			local j = `i' + 1
			mat LT_NoLLF_male[`j' , 4] 		= LT_NoLLF_male[`i' , 4] - LT_NoLLF_male[`i',5]
			mat LT_NoLLF_female[`j' , 4] 	= LT_NoLLF_female[`i' , 4] - LT_NoLLF_female[`i',5]
			mat LT_LLF_male[`j' , 4] 		= LT_LLF_male[`i' , 4] - LT_LLF_male[`i',5]
			mat LT_LLF_female[`j' , 4] 		= LT_LLF_female[`i' , 4] - LT_LLF_female[`i',5]
		} 
	}
		
	mat colnames LT_NoLLF_male 		= "Month" "XB_NoLLF_male" "Qx_NoLLF_male" "Lx_NoLLF_male" "Dx_NoLLF_male" 
	mat colnames LT_NoLLF_female 	= "Month" "XB_NoLLF_female" "Qx_NoLLF_female" "Lx_NoLLF_female" "Dx_NoLLF_female" 
	mat colnames LT_LLF_male 		= "Month" "XB_LLF_male" "Qx_LLF_male" "Lx_LLF_male" "Dx_LLF_male"
	mat colnames LT_LLF_female 		= "Month" "XB_LLF_female" "Qx_LLF_female" "Lx_LLF_female" "Dx_LLF_female"
	
	mat temp1 = LT_NoLLF_female[1..., 2...]	
	mat temp2 = LT_LLF_male[1..., 2...]	
	mat temp3 = LT_LLF_female[1..., 2...]
	
	mat colnames temp1 	= "XB_NoLLF_female" "Qx_NoLLF_female" "Lx_NoLLF_female" "Dx_NoLLF_female" 
	mat colnames temp2 	= "XB_LLF_male" "Qx_LLF_male" "Lx_LLF_male" "Dx_LLF_male"
	mat colnames temp3 	= "XB_LLF_female" "Qx_LLF_female" "Lx_LLF_female" "Dx_LLF_female"
	
	mat LifeTable_2_NoSon = LT_NoLLF_male, temp1, temp2, temp3
	 
	* Plot the CDF
	clear 
	svmat LifeTable_2_NoSon, names(col)

	* Compute the change in mean age and median interval length (quarters) 
	local mean_NoLLF_male = 0
	local dx_cum_NoLLF_male = 0 
	local flag_NoLLF_male = 1

	forval i = 1/26 { 
		local j = `i' - 1
		local int = LT_NoLLF_male[`i', 1]*LT_NoLLF_male[`i', 5]
		local dx = LT_NoLLF_male[`i', 5]
		local mean_NoLLF_male = `mean_NoLLF_male' + `int' 
		local dx_cum_NoLLF_male = `dx_cum_NoLLF_male' + `dx' 
		
		local y = LT_NoLLF_male[`i', 4]
		if `y' < .5 & `flag_NoLLF_male' == 1 { 
			local flag_NoLLF_male = 2 
			local y_1 = LT_NoLLF_male[`j', 4]
			local x = LT_NoLLF_male[`i', 1]
			local beta = `y' - `y_1'
			local alpha = `y'-(`beta'*`x')
			local median_NoLLF_male = (.5-`alpha')/`beta'
			}
	}		

	local mean_NoLLF_female = 0
	local dx_cum_NoLLF_female = 0 
	local flag_NoLLF_female = 1

	forval i = 1/26 { 
		local j = `i' - 1
		local int = LT_NoLLF_female[`i', 1]*LT_NoLLF_female[`i', 5]
		local dx = LT_NoLLF_female[`i', 5]
		local mean_NoLLF_female = `mean_NoLLF_female' + `int' 
		local dx_cum_NoLLF_female = `dx_cum_NoLLF_female' + `dx' 
		
		local y = LT_NoLLF_female[`i', 4]
		if `y' < .5 & `flag_NoLLF_female' == 1 { 
			local flag_NoLLF_female = 2 
			local y_1 = LT_NoLLF_female[`j', 4]
			local x = LT_NoLLF_female[`i', 1]
			local beta = `y' - `y_1'
			local alpha = `y'-(`beta'*`x')
			local median_NoLLF_female = (.5-`alpha')/`beta'
			}
	}		
	
	local mean_LLF_male = 0 
	local dx_cum_LLF_male = 0 
	local flag_LLF_male = 1
	forval i = 1/28 { 
		local j = `i' - 1
		local int = LT_LLF_male[`i', 1]*LT_LLF_male[`i', 5]
		local dx = LT_LLF_male[`i', 5]
		local mean_LLF_male = `mean_LLF_male' + `int' 
		local dx_cum_LLF_male = `dx_cum_LLF_male' + `dx' 
		
		local y = LT_LLF_male[`i', 4]
		if `y' < .5 & `flag_LLF_male' == 1 { 
			local flag_LLF_male = 2 
			local y_1 = LT_LLF_male[`j', 4]
			local x = LT_LLF_male[`i', 1]
			local beta = `y' - `y_1'
			local alpha = `y'-(`beta'*`x')
			local median_LLF_male = (.5-`alpha')/`beta'
			}
		} 

	local mean_LLF_female = 0 
	local dx_cum_LLF_female = 0 
	local flag_LLF_female = 1
	forval i = 1/28 { 
		local j = `i' - 1
		local int = LT_LLF_female[`i', 1]*LT_LLF_female[`i', 5]
		local dx = LT_LLF_female[`i', 5]
		local mean_LLF_female = `mean_LLF_female' + `int' 
		local dx_cum_LLF_female = `dx_cum_LLF_female' + `dx' 
		
		local y = LT_LLF_female[`i', 4]
		if `y' < .5 & `flag_LLF_female' == 1 { 
			local flag_LLF_female = 2 
			local y_1 = LT_LLF_female[`j', 4]
			local x = LT_LLF_female[`i', 1]
			local beta = `y' - `y_1'
			local alpha = `y'-(`beta'*`x')
			local median_LLF_female = (.5-`alpha')/`beta'
			}
		} 
		
	local mean_NoLLF_male 	= `mean_NoLLF_male'*`dx_cum_NoLLF_male'
	local mean_NoLLF_female = `mean_NoLLF_female'*`dx_cum_NoLLF_female'
	local mean_LLF_male 	= `mean_LLF_male'*`dx_cum_LLF_male'
	local mean_LLF_female 	= `mean_LLF_female'*`dx_cum_LLF_female'
	
	local change_in_gap = (`mean_LLF_male'-`mean_LLF_female')-(`mean_NoLLF_male'-`mean_NoLLF_female')
	
	mat BF_effects_Second_NoSon = `change_in_gap'

	restore 
	
	***************************************
	* (2.3) Second With Sons
	***************************************
			
	preserve
		
	keep if parity == 1  & year_marriage >= 1960 & birth_year <= 1979 & son == 1
	keep marriage_age cmale m_year_birth educ U5MR_5yr_ave $controls2 birth_month birth_year breastfeed_month LLFyear year uniqueid prov svydate marcmc
		
	gen 	month0 = 1 if breastfeed_month == 1 & breastfeed_month != .
	
	forval i = 1/60 { 
		gen 	month`i' = 0 if breastfeed_month > `i' & breastfeed_month != .
		replace month`i' = 1 if breastfeed_month == `i'  	
	}
		
	replace month60 = 1 if breastfeed_month > 60 & breastfeed_month != .
		
	gen start_cmc = 12*(birth_year-1900) + birth_month
		
		forval i = 1/60 { 
			cap drop cmc`i'
			gen cmc`i' = .
			
			replace cmc`i' = start_cmc + `i' if month`i' != .
			
			cap drop pd_year`i'  
			gen pd_year`i' = . 
			replace pd_year`i' = cmc`i'/12 + 1900 if month`i' != .
			
			cap drop pd_LLF`i' 
			gen pd_LLF`i' = . 
			replace pd_LLF`i' = 0 if month`i' != . 
			replace pd_LLF`i' = 1 if pd_year`i' >= LLFyear & month`i' != . 
			} 
			
	drop cmc* start_cmc
		
	reshape long pd_LLF pd_year month, i(uniqueid) j(interval_month)
	replace pd_year = floor(pd_year)
	rename month end_bf
	rename interval_month month
	stset month, failure(end_bf) id(uniqueid)
	drop if end_bf==.
			
	tab month, gen(month_)
	forval i = 1/60 { 
		gen LLF_month_`i' = pd_LLF*month_`i' 
		} 

	forval i = 1/60 { 
		gen month_`i'_male = month_`i'*cmale 
		} 		
		
	gen pd_LLF_male = pd_LLF*cmale
		
	forval i = 1/60 { 
		gen LLF_month_`i'_male = pd_LLF*month_`i'*cmale 
		} 
			
	gen birth_year2LLF = birth_year-LLFyear
		
	keep if pd_year >= 1964
	keep if birth_year2LLF >= -4	
	keep if birth_year2LLF <= 9		
		
	xi: logit 	end_bf 	pd_LLF cmale month_1-month_11 month_13-month_60 ///
						LLF_month_1-LLF_month_11 LLF_month_13-LLF_month_60 pd_LLF_male ///
						month_1_male-month_11_male month_13_male-month_60_male /// 
						LLF_month_1_male-LLF_month_11_male LLF_month_13_male-LLF_month_60_male ///
						m_year_birth $controls1 $controls2 ///
						i.prov i.pd_year if pd_year!=1967 & pd_year != 1968, cluster(prov) ltolerance(1e-1) iter(30)

	* REFERENCE GROUP: MALES PRE LLF
	replace pd_LLF = 0
	replace cmale = 1 

	forval i = 1/60 { 			
		replace LLF_month_`i' = 0
		} 
			
	replace pd_LLF_male = 0 
		
	forval i = 1/60 { 			
		replace month_`i'_male = month_`i'*cmale
		}

	forval i = 1/60 { 			
		replace LLF_month_`i'_male = 0
		}
		
	cap drop xb
	predict xb, xb
		
	cap drop qx
	predict qx
		
		forval q = 1/60 { 
			
			sum xb if month == `q' 
			local xb = r(mean)
				
			sum qx if month == `q' 
			local qx = r(mean)
				
			if `q' == 1 {
				mat LT_NoLLF_male = `q', `xb', `qx' 
				} 
			else if `q' > 1 { 
				mat temp = `q', `xb', `qx' 
				mat LT_NoLLF_male = LT_NoLLF_male \ temp
				}
			}
				
	* FEMALES PRE LLF
	replace cmale = 0 
		
	forval i = 1/60 { 			
		replace month_`i'_male = 0
		}
		
	cap drop xb
	predict xb, xb
		
	cap drop qx
	predict qx
		
		forval q = 1/60 { 
			
			sum xb if month == `q' 
			local xb = r(mean)
				
			sum qx if month == `q' 
			local qx = r(mean)
				
			if `q' == 1 {
				mat LT_NoLLF_female = `q', `xb', `qx' 
				} 
			else if `q' > 1 { 
				mat temp = `q', `xb', `qx' 
				mat LT_NoLLF_female = LT_NoLLF_female \ temp
				}
			}

		
	* MALES POST LLF
	replace pd_LLF = 1
	replace cmale = 1
		
	forval i = 1/60 { 			
		replace LLF_month_`i' = month_`i'*pd_LLF 
		} 
			
	replace pd_LLF_male = 1 
	
	forval i = 1/60 { 			
		replace month_`i'_male = month_`i'*cmale
		}
			
	forval i = 1/60 { 			
		replace LLF_month_`i'_male = month_`i'*cmale*pd_LLF
		}
			
	cap drop xb
	predict xb, xb
		
	cap drop qx
	predict qx
		
			forval q = 1/60 { 
			
				sum xb if month == `q' 
				local xb = r(mean)
				
				sum qx if month == `q' 
				local qx = r(mean)
				
				if `q' == 1 {
					mat LT_LLF_male = `q', `xb', `qx' 
					} 
				else if `q' > 1 { 
					mat temp = `q', `xb', `qx' 
					mat LT_LLF_male = LT_LLF_male \ temp
					}
				}

	* FEMALES POST LLF
	replace cmale = 0
	replace pd_LLF_male = 0
		
	forval i = 1/60 { 			
		replace month_`i'_male = 0
		}
			
	forval i = 1/60 { 			
		replace LLF_month_`i'_male = 0
		}
			
	cap drop xb
	predict xb, xb
		
	cap drop qx
	predict qx
		
			forval q = 1/60 { 
			
				sum xb if month == `q' 
				local xb = r(mean)
				
				sum qx if month == `q' 
				local qx = r(mean)
				
				if `q' == 1 {
					mat LT_LLF_female = `q', `xb', `qx' 
					} 
				else if `q' > 1 { 
					mat temp = `q', `xb', `qx' 
					mat LT_LLF_female = LT_LLF_female \ temp
					}
				}
				
	* Compute Lx and Dx for each month, males and females, with and without LLF
	mat temp = J(60, 1, .) 
	mat LT_NoLLF_male 	= LT_NoLLF_male, temp, temp
	mat LT_NoLLF_female = LT_NoLLF_female, temp, temp
	mat LT_LLF_male 	= LT_LLF_male, temp, temp
	mat LT_LLF_female 	= LT_LLF_female, temp, temp
		
	mat LT_NoLLF_male[1,4] 		= 1
	mat LT_NoLLF_female[1,4] 	= 1
	mat LT_LLF_male[1,4] 		= 1
	mat LT_LLF_female[1,4] 		= 1
	
	forval i = 1/60 {
		mat LT_NoLLF_male[`i',5] 	= LT_NoLLF_male[`i',3]*LT_NoLLF_male[`i',4]
		mat LT_NoLLF_female[`i',5] 	= LT_NoLLF_female[`i',3]*LT_NoLLF_female[`i',4]
		mat LT_LLF_male[`i',5] 		= LT_LLF_male[`i',3]*LT_LLF_male[`i',4]
		mat LT_LLF_female[`i',5] 	= LT_LLF_female[`i',3]*LT_LLF_female[`i',4]
		
		if `i' < 60 { 
			local j = `i' + 1
			mat LT_NoLLF_male[`j' , 4] 		= LT_NoLLF_male[`i' , 4] - LT_NoLLF_male[`i',5]
			mat LT_NoLLF_female[`j' , 4] 	= LT_NoLLF_female[`i' , 4] - LT_NoLLF_female[`i',5]
			mat LT_LLF_male[`j' , 4] 		= LT_LLF_male[`i' , 4] - LT_LLF_male[`i',5]
			mat LT_LLF_female[`j' , 4] 		= LT_LLF_female[`i' , 4] - LT_LLF_female[`i',5]
		} 
	}
		
	mat colnames LT_NoLLF_male 		= "Month" "XB_NoLLF_male" "Qx_NoLLF_male" "Lx_NoLLF_male" "Dx_NoLLF_male" 
	mat colnames LT_NoLLF_female 	= "Month" "XB_NoLLF_female" "Qx_NoLLF_female" "Lx_NoLLF_female" "Dx_NoLLF_female" 
	mat colnames LT_LLF_male 		= "Month" "XB_LLF_male" "Qx_LLF_male" "Lx_LLF_male" "Dx_LLF_male"
	mat colnames LT_LLF_female 		= "Month" "XB_LLF_female" "Qx_LLF_female" "Lx_LLF_female" "Dx_LLF_female"
	
	mat temp1 = LT_NoLLF_female[1..., 2...]	
	mat temp2 = LT_LLF_male[1..., 2...]	
	mat temp3 = LT_LLF_female[1..., 2...]
	
	mat colnames temp1 	= "XB_NoLLF_female" "Qx_NoLLF_female" "Lx_NoLLF_female" "Dx_NoLLF_female" 
	mat colnames temp2 	= "XB_LLF_male" "Qx_LLF_male" "Lx_LLF_male" "Dx_LLF_male"
	mat colnames temp3 	= "XB_LLF_female" "Qx_LLF_female" "Lx_LLF_female" "Dx_LLF_female"
	
	mat LifeTable_2_Son = LT_NoLLF_male, temp1, temp2, temp3
	 
	clear 
	svmat LifeTable_2_Son, names(col)
	
	* Compute the change in mean age and median interval length (quarters) 
	local mean_NoLLF_male = 0
	local dx_cum_NoLLF_male = 0 
	local flag_NoLLF_male = 1

	forval i = 1/26 { 
		local j = `i' - 1
		local int = LT_NoLLF_male[`i', 1]*LT_NoLLF_male[`i', 5]
		local dx = LT_NoLLF_male[`i', 5]
		local mean_NoLLF_male = `mean_NoLLF_male' + `int' 
		local dx_cum_NoLLF_male = `dx_cum_NoLLF_male' + `dx' 
		
		local y = LT_NoLLF_male[`i', 4]
		if `y' < .5 & `flag_NoLLF_male' == 1 { 
			local flag_NoLLF_male = 2 
			local y_1 = LT_NoLLF_male[`j', 4]
			local x = LT_NoLLF_male[`i', 1]
			local beta = `y' - `y_1'
			local alpha = `y'-(`beta'*`x')
			local median_NoLLF_male = (.5-`alpha')/`beta'
			}
	}		

	local mean_NoLLF_female = 0
	local dx_cum_NoLLF_female = 0 
	local flag_NoLLF_female = 1

	forval i = 1/26 { 
		local j = `i' - 1
		local int = LT_NoLLF_female[`i', 1]*LT_NoLLF_female[`i', 5]
		local dx = LT_NoLLF_female[`i', 5]
		local mean_NoLLF_female = `mean_NoLLF_female' + `int' 
		local dx_cum_NoLLF_female = `dx_cum_NoLLF_female' + `dx' 
		
		local y = LT_NoLLF_female[`i', 4]
		if `y' < .5 & `flag_NoLLF_female' == 1 { 
			local flag_NoLLF_female = 2 
			local y_1 = LT_NoLLF_female[`j', 4]
			local x = LT_NoLLF_female[`i', 1]
			local beta = `y' - `y_1'
			local alpha = `y'-(`beta'*`x')
			local median_NoLLF_female = (.5-`alpha')/`beta'
			}
	}		
	
	local mean_LLF_male = 0 
	local dx_cum_LLF_male = 0 
	local flag_LLF_male = 1
	forval i = 1/28 { 
		local j = `i' - 1
		local int = LT_LLF_male[`i', 1]*LT_LLF_male[`i', 5]
		local dx = LT_LLF_male[`i', 5]
		local mean_LLF_male = `mean_LLF_male' + `int' 
		local dx_cum_LLF_male = `dx_cum_LLF_male' + `dx' 
		
		local y = LT_LLF_male[`i', 4]
		if `y' < .5 & `flag_LLF_male' == 1 { 
			local flag_LLF_male = 2 
			local y_1 = LT_LLF_male[`j', 4]
			local x = LT_LLF_male[`i', 1]
			local beta = `y' - `y_1'
			local alpha = `y'-(`beta'*`x')
			local median_LLF_male = (.5-`alpha')/`beta'
			}
		} 

	local mean_LLF_female = 0 
	local dx_cum_LLF_female = 0 
	local flag_LLF_female = 1
	forval i = 1/28 { 
		local j = `i' - 1
		local int = LT_LLF_female[`i', 1]*LT_LLF_female[`i', 5]
		local dx = LT_LLF_female[`i', 5]
		local mean_LLF_female = `mean_LLF_female' + `int' 
		local dx_cum_LLF_female = `dx_cum_LLF_female' + `dx' 
		
		local y = LT_LLF_female[`i', 4]
		if `y' < .5 & `flag_LLF_female' == 1 { 
			local flag_LLF_female = 2 
			local y_1 = LT_LLF_female[`j', 4]
			local x = LT_LLF_female[`i', 1]
			local beta = `y' - `y_1'
			local alpha = `y'-(`beta'*`x')
			local median_LLF_female = (.5-`alpha')/`beta'
			}
		} 
			
	local mean_NoLLF_male 	= `mean_NoLLF_male'*`dx_cum_NoLLF_male'
	local mean_NoLLF_female = `mean_NoLLF_female'*`dx_cum_NoLLF_female'
	local mean_LLF_male 	= `mean_LLF_male'*`dx_cum_LLF_male'
	local mean_LLF_female 	= `mean_LLF_female'*`dx_cum_LLF_female'
	
	local change_in_gap = (`mean_LLF_male'-`mean_LLF_female')-(`mean_NoLLF_male'-`mean_NoLLF_female')
	
	mat BF_effects_Second_Son = `change_in_gap'
	
	restore 
	
	*******************************************
	* (2.4) Third And Higher, No Previous Sons
	*******************************************
			
	preserve
		
	keep if parity >= 2 & year_marriage >= 1960 & birth_year <= 1979 & son == 0
	keep marriage_age parity cmale m_year_birth educ U5MR_5yr_ave $controls2 birth_month birth_year breastfeed_month LLFyear year uniqueid prov svydate marcmc
		
	gen 	month0 = 1 if breastfeed_month == 1 & breastfeed_month != .
		
		forval i = 1/60 { 
			gen 	month`i' = 0 if breastfeed_month > `i' & breastfeed_month != .
			replace month`i' = 1 if breastfeed_month == `i'  	
		}
			
	replace month60 = 1 if breastfeed_month > 60 & breastfeed_month != .
		
	gen start_cmc = 12*(birth_year-1900) + birth_month
		
		forval i = 1/60 { 
			cap drop cmc`i'
			gen cmc`i' = .
			
			replace cmc`i' = start_cmc + `i' if month`i' != .
			
			cap drop pd_year`i'  
			gen pd_year`i' = . 
			replace pd_year`i' = cmc`i'/12 + 1900 if month`i' != .
			
			cap drop pd_LLF`i' 
			gen pd_LLF`i' = . 
			replace pd_LLF`i' = 0 if month`i' != . 
			replace pd_LLF`i' = 1 if pd_year`i' >= LLFyear & month`i' != . 
			} 
			
	drop cmc* start_cmc
		
	reshape long pd_LLF pd_year month, i(uniqueid parity birth_year birth_month) j(interval_month)
	replace pd_year = floor(pd_year)
	rename month end_bf
	rename interval_month month
	stset month, failure(end_bf) id(uniqueid)
	drop if end_bf==.
			
	tab month, gen(month_)
		forval i = 1/60 { 
			gen LLF_month_`i' = pd_LLF*month_`i' 
			} 

		forval i = 1/60 { 
			gen month_`i'_male = month_`i'*cmale 
			} 		
		
		gen pd_LLF_male = pd_LLF*cmale
		
		forval i = 1/60 { 
			gen LLF_month_`i'_male = pd_LLF*month_`i'*cmale 
			} 
			
	gen birth_year2LLF = birth_year-LLFyear
		
	keep if pd_year >= 1964
	keep if birth_year2LLF >= -4	
	keep if birth_year2LLF <= 9		
				
	xi: logit 	end_bf 	pd_LLF cmale month_1-month_11 month_13-month_60 ///
						LLF_month_1-LLF_month_11 LLF_month_13-LLF_month_60 pd_LLF_male ///
						month_1_male-month_11_male month_13_male-month_60_male /// 
						LLF_month_1_male-LLF_month_11_male LLF_month_13_male-LLF_month_60_male ///
						m_year_birth $controls1 $controls2 ///
						i.prov i.pd_year if pd_year!=1967 & pd_year != 1968, cluster(prov) ltolerance(1e-1) iter(30)

	* REFERENCE GROUP: MALES PRE LLF
	replace pd_LLF = 0
	replace cmale = 1 

	forval i = 1/60 { 			
		replace LLF_month_`i' = 0
		} 
			
	replace pd_LLF_male = 0 
		
	forval i = 1/60 { 			
		replace month_`i'_male = month_`i'*cmale
		}

	forval i = 1/60 { 			
		replace LLF_month_`i'_male = 0
		}
		
	cap drop xb
	predict xb, xb
		
	cap drop qx
	predict qx
		
			forval q = 1/60 { 
			
				sum xb if month == `q' 
				local xb = r(mean)
				
				sum qx if month == `q' 
				local qx = r(mean)
				
				if `q' == 1 {
					mat LT_NoLLF_male = `q', `xb', `qx' 
					} 
				else if `q' > 1 { 
					mat temp = `q', `xb', `qx' 
					mat LT_NoLLF_male = LT_NoLLF_male \ temp
					}
				}
				
	* FEMALES PRE LLF
	replace cmale = 0 
		
	forval i = 1/60 { 			
		replace month_`i'_male = 0
		}
		
	cap drop xb
	predict xb, xb
		
	cap drop qx
	predict qx
		
		forval q = 1/60 { 
			
			sum xb if month == `q' 
			local xb = r(mean)
				
			sum qx if month == `q' 
			local qx = r(mean)
				
			if `q' == 1 {
				mat LT_NoLLF_female = `q', `xb', `qx' 
				} 
			else if `q' > 1 { 
				mat temp = `q', `xb', `qx' 
				mat LT_NoLLF_female = LT_NoLLF_female \ temp
				}
			}

	* MALES POST LLF
	replace pd_LLF = 1	
	replace cmale = 1
		
	forval i = 1/60 { 			
		replace LLF_month_`i' = month_`i'*pd_LLF 
		} 
			
	replace pd_LLF_male = 1 
		
	forval i = 1/60 { 			
		replace month_`i'_male = month_`i'*cmale
		}
			
	forval i = 1/60 { 			
		replace LLF_month_`i'_male = month_`i'*cmale*pd_LLF
		}
			
	cap drop xb
	predict xb, xb
		
	cap drop qx
	predict qx
		
			forval q = 1/60 { 
			
				sum xb if month == `q' 
				local xb = r(mean)
				
				sum qx if month == `q' 
				local qx = r(mean)
				
				if `q' == 1 {
					mat LT_LLF_male = `q', `xb', `qx' 
					} 
				else if `q' > 1 { 
					mat temp = `q', `xb', `qx' 
					mat LT_LLF_male = LT_LLF_male \ temp
					}
				}

		
	* FEMALES POST LLF
	replace cmale = 0	
	replace pd_LLF_male = 0
		
	forval i = 1/60 { 			
		replace month_`i'_male = 0
		}
			
	forval i = 1/60 { 			
		replace LLF_month_`i'_male = 0
		}
			
	cap drop xb
	predict xb, xb
		
	cap drop qx
	predict qx
		
			forval q = 1/60 { 
			
				sum xb if month == `q' 
				local xb = r(mean)
				
				sum qx if month == `q' 
				local qx = r(mean)
				
				if `q' == 1 {
					mat LT_LLF_female = `q', `xb', `qx' 
					} 
				else if `q' > 1 { 
					mat temp = `q', `xb', `qx' 
					mat LT_LLF_female = LT_LLF_female \ temp
					}
				}
				
	mat temp = J(60, 1, .) 
	mat LT_NoLLF_male 	= LT_NoLLF_male, temp, temp
	mat LT_NoLLF_female = LT_NoLLF_female, temp, temp
	mat LT_LLF_male 	= LT_LLF_male, temp, temp
	mat LT_LLF_female 	= LT_LLF_female, temp, temp
		
	mat LT_NoLLF_male[1,4] 		= 1
	mat LT_NoLLF_female[1,4] 	= 1
	mat LT_LLF_male[1,4] 		= 1
	mat LT_LLF_female[1,4] 		= 1
	
	forval i = 1/60 {
		mat LT_NoLLF_male[`i',5] 	= LT_NoLLF_male[`i',3]*LT_NoLLF_male[`i',4]
		mat LT_NoLLF_female[`i',5] 	= LT_NoLLF_female[`i',3]*LT_NoLLF_female[`i',4]
		mat LT_LLF_male[`i',5] 		= LT_LLF_male[`i',3]*LT_LLF_male[`i',4]
		mat LT_LLF_female[`i',5] 	= LT_LLF_female[`i',3]*LT_LLF_female[`i',4]
		
		if `i' < 60 { 
			local j = `i' + 1
			mat LT_NoLLF_male[`j' , 4] 		= LT_NoLLF_male[`i' , 4] - LT_NoLLF_male[`i',5]
			mat LT_NoLLF_female[`j' , 4] 	= LT_NoLLF_female[`i' , 4] - LT_NoLLF_female[`i',5]
			mat LT_LLF_male[`j' , 4] 		= LT_LLF_male[`i' , 4] - LT_LLF_male[`i',5]
			mat LT_LLF_female[`j' , 4] 		= LT_LLF_female[`i' , 4] - LT_LLF_female[`i',5]
		} 
	}
		
	mat colnames LT_NoLLF_male 		= "Month" "XB_NoLLF_male" "Qx_NoLLF_male" "Lx_NoLLF_male" "Dx_NoLLF_male" 
	mat colnames LT_NoLLF_female 	= "Month" "XB_NoLLF_female" "Qx_NoLLF_female" "Lx_NoLLF_female" "Dx_NoLLF_female" 
	mat colnames LT_LLF_male 		= "Month" "XB_LLF_male" "Qx_LLF_male" "Lx_LLF_male" "Dx_LLF_male"
	mat colnames LT_LLF_female 		= "Month" "XB_LLF_female" "Qx_LLF_female" "Lx_LLF_female" "Dx_LLF_female"
	
	mat temp1 = LT_NoLLF_female[1..., 2...]	
	mat temp2 = LT_LLF_male[1..., 2...]	
	mat temp3 = LT_LLF_female[1..., 2...]
	
	mat colnames temp1 	= "XB_NoLLF_female" "Qx_NoLLF_female" "Lx_NoLLF_female" "Dx_NoLLF_female" 
	mat colnames temp2 	= "XB_LLF_male" "Qx_LLF_male" "Lx_LLF_male" "Dx_LLF_male"
	mat colnames temp3 	= "XB_LLF_female" "Qx_LLF_female" "Lx_LLF_female" "Dx_LLF_female"
	
	mat LifeTable_3_NoSon = LT_NoLLF_male, temp1, temp2, temp3
	 	
	clear 
	svmat LifeTable_3_NoSon, names(col)
	
	* Compute the change in mean age and median interval length (quarters) 
	local mean_NoLLF_male = 0
	local dx_cum_NoLLF_male = 0 
	local flag_NoLLF_male = 1

	forval i = 1/26 { 
		local j = `i' - 1
		local int = LT_NoLLF_male[`i', 1]*LT_NoLLF_male[`i', 5]
		local dx = LT_NoLLF_male[`i', 5]
		local mean_NoLLF_male = `mean_NoLLF_male' + `int' 
		local dx_cum_NoLLF_male = `dx_cum_NoLLF_male' + `dx' 
		
		local y = LT_NoLLF_male[`i', 4]
		if `y' < .5 & `flag_NoLLF_male' == 1 { 
			local flag_NoLLF_male = 2 
			local y_1 = LT_NoLLF_male[`j', 4]
			local x = LT_NoLLF_male[`i', 1]
			local beta = `y' - `y_1'
			local alpha = `y'-(`beta'*`x')
			local median_NoLLF_male = (.5-`alpha')/`beta'
			}
	}		

	local mean_NoLLF_female = 0
	local dx_cum_NoLLF_female = 0 
	local flag_NoLLF_female = 1

	forval i = 1/26 { 
		local j = `i' - 1
		local int = LT_NoLLF_female[`i', 1]*LT_NoLLF_female[`i', 5]
		local dx = LT_NoLLF_female[`i', 5]
		local mean_NoLLF_female = `mean_NoLLF_female' + `int' 
		local dx_cum_NoLLF_female = `dx_cum_NoLLF_female' + `dx' 
		
		local y = LT_NoLLF_female[`i', 4]
		if `y' < .5 & `flag_NoLLF_female' == 1 { 
			local flag_NoLLF_female = 2 
			local y_1 = LT_NoLLF_female[`j', 4]
			local x = LT_NoLLF_female[`i', 1]
			local beta = `y' - `y_1'
			local alpha = `y'-(`beta'*`x')
			local median_NoLLF_female = (.5-`alpha')/`beta'
			}
	}		
	
	local mean_LLF_male = 0 
	local dx_cum_LLF_male = 0 
	local flag_LLF_male = 1
	forval i = 1/28 { 
		local j = `i' - 1
		local int = LT_LLF_male[`i', 1]*LT_LLF_male[`i', 5]
		local dx = LT_LLF_male[`i', 5]
		local mean_LLF_male = `mean_LLF_male' + `int' 
		local dx_cum_LLF_male = `dx_cum_LLF_male' + `dx' 
		
		local y = LT_LLF_male[`i', 4]
		if `y' < .5 & `flag_LLF_male' == 1 { 
			local flag_LLF_male = 2 
			local y_1 = LT_LLF_male[`j', 4]
			local x = LT_LLF_male[`i', 1]
			local beta = `y' - `y_1'
			local alpha = `y'-(`beta'*`x')
			local median_LLF_male = (.5-`alpha')/`beta'
			}
		} 

	local mean_LLF_female = 0 
	local dx_cum_LLF_female = 0 
	local flag_LLF_female = 1
	forval i = 1/28 { 
		local j = `i' - 1
		local int = LT_LLF_female[`i', 1]*LT_LLF_female[`i', 5]
		local dx = LT_LLF_female[`i', 5]
		local mean_LLF_female = `mean_LLF_female' + `int' 
		local dx_cum_LLF_female = `dx_cum_LLF_female' + `dx' 
		
		local y = LT_LLF_female[`i', 4]
		if `y' < .5 & `flag_LLF_female' == 1 { 
			local flag_LLF_female = 2 
			local y_1 = LT_LLF_female[`j', 4]
			local x = LT_LLF_female[`i', 1]
			local beta = `y' - `y_1'
			local alpha = `y'-(`beta'*`x')
			local median_LLF_female = (.5-`alpha')/`beta'
			}
		} 
		
	local mean_NoLLF_male 	= `mean_NoLLF_male'*`dx_cum_NoLLF_male'
	local mean_NoLLF_female = `mean_NoLLF_female'*`dx_cum_NoLLF_female'
	local mean_LLF_male 	= `mean_LLF_male'*`dx_cum_LLF_male'
	local mean_LLF_female 	= `mean_LLF_female'*`dx_cum_LLF_female'
	
	local change_in_gap = (`mean_LLF_male'-`mean_LLF_female')-(`mean_NoLLF_male'-`mean_NoLLF_female')
	
	mat BF_effects_Third_NoSon = `change_in_gap'

	restore

	***************************************
	* (2.5) Third And Higher With Sons
	***************************************
			
	preserve
		
	keep if parity >= 2  & year_marriage >= 1960 & birth_year <= 1979 & son == 1
	keep marriage_age parity cmale m_year_birth educ U5MR_5yr_ave $controls2 birth_month birth_year breastfeed_month LLFyear year uniqueid prov svydate marcmc
	
	gen 	month0 = 1 if breastfeed_month == 1 & breastfeed_month != .
		
		forval i = 1/60 { 
			gen 	month`i' = 0 if breastfeed_month > `i' & breastfeed_month != .
			replace month`i' = 1 if breastfeed_month == `i'  	
		}
		
	replace month60 = 1 if breastfeed_month > 60 & breastfeed_month != .
		
	gen start_cmc = 12*(birth_year-1900) + birth_month
		
		forval i = 1/60 { 
			cap drop cmc`i'
			gen cmc`i' = .
			
			replace cmc`i' = start_cmc + `i' if month`i' != .
			
			cap drop pd_year`i'  
			gen pd_year`i' = . 
			replace pd_year`i' = cmc`i'/12 + 1900 if month`i' != .
			
			cap drop pd_LLF`i' 
			gen pd_LLF`i' = . 
			replace pd_LLF`i' = 0 if month`i' != . 
			replace pd_LLF`i' = 1 if pd_year`i' >= LLFyear & month`i' != . 
			} 
			
	drop cmc* start_cmc
		
	reshape long pd_LLF pd_year month, i(uniqueid parity birth_year birth_month) j(interval_month)
	replace pd_year = floor(pd_year)
	rename month end_bf
	rename interval_month month
	stset month, failure(end_bf) id(uniqueid)
	drop if end_bf==.
			
	tab month, gen(month_)
	forval i = 1/60 { 
		gen LLF_month_`i' = pd_LLF*month_`i' 
		} 

	forval i = 1/60 { 
		gen month_`i'_male = month_`i'*cmale 
		} 		
		
	gen pd_LLF_male = pd_LLF*cmale
		
	forval i = 1/60 { 
		gen LLF_month_`i'_male = pd_LLF*month_`i'*cmale 
		} 
			
	gen birth_year2LLF = birth_year-LLFyear
		
	keep if pd_year >= 1964
	keep if birth_year2LLF >= -4	
	keep if birth_year2LLF <= 9		
		
	xi: logit 	end_bf 	pd_LLF cmale month_1-month_11 month_13-month_60 ///
						LLF_month_1-LLF_month_11 LLF_month_13-LLF_month_60 pd_LLF_male ///
						month_1_male-month_11_male month_13_male-month_60_male /// 
						LLF_month_1_male-LLF_month_11_male LLF_month_13_male-LLF_month_60_male ///
						m_year_birth $controls1 $controls2 ///
						i.prov i.pd_year if pd_year!=1967 & pd_year != 1968, cluster(prov) ltolerance(1e-1) iter(30)

	* REFERENCE GROUP: MALES PRE LLF
	replace pd_LLF = 0
	replace cmale = 1 

		forval i = 1/60 { 			
			replace LLF_month_`i' = 0
			} 
			
	replace pd_LLF_male = 0 
		
		forval i = 1/60 { 			
			replace month_`i'_male = month_`i'*cmale
			}

		forval i = 1/60 { 			
			replace LLF_month_`i'_male = 0
			}
		
	cap drop xb
	predict xb, xb
		
	cap drop qx
	predict qx
		
		forval q = 1/60 { 
			
			sum xb if month == `q' 
			local xb = r(mean)
				
			sum qx if month == `q' 
			local qx = r(mean)
				
				if `q' == 1 {
					mat LT_NoLLF_male = `q', `xb', `qx' 
					} 
				else if `q' > 1 { 
					mat temp = `q', `xb', `qx' 
					mat LT_NoLLF_male = LT_NoLLF_male \ temp
					}
				}
				
	* FEMALES PRE LLF
	replace cmale = 0 
		
		forval i = 1/60 { 			
			replace month_`i'_male = 0
			}
		
	cap drop xb
	predict xb, xb
		
	cap drop qx
	predict qx
		
			forval q = 1/60 { 
			
				sum xb if month == `q' 
				local xb = r(mean)
				
				sum qx if month == `q' 
				local qx = r(mean)
				
				if `q' == 1 {
					mat LT_NoLLF_female = `q', `xb', `qx' 
					} 
				else if `q' > 1 { 
					mat temp = `q', `xb', `qx' 
					mat LT_NoLLF_female = LT_NoLLF_female \ temp
					}
				}
		
	* MALES POST LLF
	replace pd_LLF = 1
	replace cmale = 1
		
		forval i = 1/60 { 			
			replace LLF_month_`i' = month_`i'*pd_LLF 
			} 
			
	replace pd_LLF_male = 1 
		
		forval i = 1/60 { 			
			replace month_`i'_male = month_`i'*cmale
			}
			
		forval i = 1/60 { 			
			replace LLF_month_`i'_male = month_`i'*cmale*pd_LLF
			}
			
	cap drop xb
	predict xb, xb
		
	cap drop qx
	predict qx
		
			forval q = 1/60 { 
			
				sum xb if month == `q' 
				local xb = r(mean)
				
				sum qx if month == `q' 
				local qx = r(mean)
				
				if `q' == 1 {
					mat LT_LLF_male = `q', `xb', `qx' 
					} 
				else if `q' > 1 { 
					mat temp = `q', `xb', `qx' 
					mat LT_LLF_male = LT_LLF_male \ temp
					}
				}

		
	* FEMALES POST LLF
	replace cmale = 0	
	replace pd_LLF_male = 0
		
		forval i = 1/60 { 			
			replace month_`i'_male = 0
			}
			
		forval i = 1/60 { 			
			replace LLF_month_`i'_male = 0
			}
			
	cap drop xb
	predict xb, xb
		
	cap drop qx
	predict qx
		
			forval q = 1/60 { 
			
				sum xb if month == `q' 
				local xb = r(mean)
				
				sum qx if month == `q' 
				local qx = r(mean)
				
				if `q' == 1 {
					mat LT_LLF_female = `q', `xb', `qx' 
					} 
				else if `q' > 1 { 
					mat temp = `q', `xb', `qx' 
					mat LT_LLF_female = LT_LLF_female \ temp
					}
				}
				
	* Compute Lx and Dx for each month, males and females, with and without LLF
	mat temp = J(60, 1, .) 
	mat LT_NoLLF_male 	= LT_NoLLF_male, temp, temp
	mat LT_NoLLF_female = LT_NoLLF_female, temp, temp
	mat LT_LLF_male 	= LT_LLF_male, temp, temp
	mat LT_LLF_female 	= LT_LLF_female, temp, temp
		
	mat LT_NoLLF_male[1,4] 		= 1
	mat LT_NoLLF_female[1,4] 	= 1
	mat LT_LLF_male[1,4] 		= 1
	mat LT_LLF_female[1,4] 		= 1
	
	forval i = 1/60 {
		mat LT_NoLLF_male[`i',5] 	= LT_NoLLF_male[`i',3]*LT_NoLLF_male[`i',4]
		mat LT_NoLLF_female[`i',5] 	= LT_NoLLF_female[`i',3]*LT_NoLLF_female[`i',4]
		mat LT_LLF_male[`i',5] 		= LT_LLF_male[`i',3]*LT_LLF_male[`i',4]
		mat LT_LLF_female[`i',5] 	= LT_LLF_female[`i',3]*LT_LLF_female[`i',4]
		
		if `i' < 60 { 
			local j = `i' + 1
			mat LT_NoLLF_male[`j' , 4] 		= LT_NoLLF_male[`i' , 4] - LT_NoLLF_male[`i',5]
			mat LT_NoLLF_female[`j' , 4] 	= LT_NoLLF_female[`i' , 4] - LT_NoLLF_female[`i',5]
			mat LT_LLF_male[`j' , 4] 		= LT_LLF_male[`i' , 4] - LT_LLF_male[`i',5]
			mat LT_LLF_female[`j' , 4] 		= LT_LLF_female[`i' , 4] - LT_LLF_female[`i',5]
		} 
	}
		
	mat colnames LT_NoLLF_male 		= "Month" "XB_NoLLF_male" "Qx_NoLLF_male" "Lx_NoLLF_male" "Dx_NoLLF_male" 
	mat colnames LT_NoLLF_female 	= "Month" "XB_NoLLF_female" "Qx_NoLLF_female" "Lx_NoLLF_female" "Dx_NoLLF_female" 
	mat colnames LT_LLF_male 		= "Month" "XB_LLF_male" "Qx_LLF_male" "Lx_LLF_male" "Dx_LLF_male"
	mat colnames LT_LLF_female 		= "Month" "XB_LLF_female" "Qx_LLF_female" "Lx_LLF_female" "Dx_LLF_female"
	
	mat temp1 = LT_NoLLF_female[1..., 2...]	
	mat temp2 = LT_LLF_male[1..., 2...]	
	mat temp3 = LT_LLF_female[1..., 2...]
	
	mat colnames temp1 	= "XB_NoLLF_female" "Qx_NoLLF_female" "Lx_NoLLF_female" "Dx_NoLLF_female" 
	mat colnames temp2 	= "XB_LLF_male" "Qx_LLF_male" "Lx_LLF_male" "Dx_LLF_male"
	mat colnames temp3 	= "XB_LLF_female" "Qx_LLF_female" "Lx_LLF_female" "Dx_LLF_female"
	
	mat LifeTable_3_Son = LT_NoLLF_male, temp1, temp2, temp3
	 
	clear 
	svmat LifeTable_3_Son, names(col)

	* Compute the change in mean age and median interval length (quarters) 
	local mean_NoLLF_male = 0
	local dx_cum_NoLLF_male = 0 
	local flag_NoLLF_male = 1

	forval i = 1/26 { 
		local j = `i' - 1
		local int = LT_NoLLF_male[`i', 1]*LT_NoLLF_male[`i', 5]
		local dx = LT_NoLLF_male[`i', 5]
		local mean_NoLLF_male = `mean_NoLLF_male' + `int' 
		local dx_cum_NoLLF_male = `dx_cum_NoLLF_male' + `dx' 
		
		local y = LT_NoLLF_male[`i', 4]
		if `y' < .5 & `flag_NoLLF_male' == 1 { 
			local flag_NoLLF_male = 2 
			local y_1 = LT_NoLLF_male[`j', 4]
			local x = LT_NoLLF_male[`i', 1]
			local beta = `y' - `y_1'
			local alpha = `y'-(`beta'*`x')
			local median_NoLLF_male = (.5-`alpha')/`beta'
			}
	}		

	local mean_NoLLF_female = 0
	local dx_cum_NoLLF_female = 0 
	local flag_NoLLF_female = 1

	forval i = 1/26 { 
		local j = `i' - 1
		local int = LT_NoLLF_female[`i', 1]*LT_NoLLF_female[`i', 5]
		local dx = LT_NoLLF_female[`i', 5]
		local mean_NoLLF_female = `mean_NoLLF_female' + `int' 
		local dx_cum_NoLLF_female = `dx_cum_NoLLF_female' + `dx' 
		
		local y = LT_NoLLF_female[`i', 4]
		if `y' < .5 & `flag_NoLLF_female' == 1 { 
			local flag_NoLLF_female = 2 
			local y_1 = LT_NoLLF_female[`j', 4]
			local x = LT_NoLLF_female[`i', 1]
			local beta = `y' - `y_1'
			local alpha = `y'-(`beta'*`x')
			local median_NoLLF_female = (.5-`alpha')/`beta'
			}
	}		
	
	local mean_LLF_male = 0 
	local dx_cum_LLF_male = 0 
	local flag_LLF_male = 1
	forval i = 1/28 { 
		local j = `i' - 1
		local int = LT_LLF_male[`i', 1]*LT_LLF_male[`i', 5]
		local dx = LT_LLF_male[`i', 5]
		local mean_LLF_male = `mean_LLF_male' + `int' 
		local dx_cum_LLF_male = `dx_cum_LLF_male' + `dx' 
		
		local y = LT_LLF_male[`i', 4]
		if `y' < .5 & `flag_LLF_male' == 1 { 
			local flag_LLF_male = 2 
			local y_1 = LT_LLF_male[`j', 4]
			local x = LT_LLF_male[`i', 1]
			local beta = `y' - `y_1'
			local alpha = `y'-(`beta'*`x')
			local median_LLF_male = (.5-`alpha')/`beta'
			}
		} 

	local mean_LLF_female = 0 
	local dx_cum_LLF_female = 0 
	local flag_LLF_female = 1
	forval i = 1/28 { 
		local j = `i' - 1
		local int = LT_LLF_female[`i', 1]*LT_LLF_female[`i', 5]
		local dx = LT_LLF_female[`i', 5]
		local mean_LLF_female = `mean_LLF_female' + `int' 
		local dx_cum_LLF_female = `dx_cum_LLF_female' + `dx' 
		
		local y = LT_LLF_female[`i', 4]
		if `y' < .5 & `flag_LLF_female' == 1 { 
			local flag_LLF_female = 2 
			local y_1 = LT_LLF_female[`j', 4]
			local x = LT_LLF_female[`i', 1]
			local beta = `y' - `y_1'
			local alpha = `y'-(`beta'*`x')
			local median_LLF_female = (.5-`alpha')/`beta'
			}
		} 
		
	local mean_NoLLF_male 	= `mean_NoLLF_male'*`dx_cum_NoLLF_male'
	local mean_NoLLF_female = `mean_NoLLF_female'*`dx_cum_NoLLF_female'
	local mean_LLF_male 	= `mean_LLF_male'*`dx_cum_LLF_male'
	local mean_LLF_female 	= `mean_LLF_female'*`dx_cum_LLF_female'
	
	
	local change_in_gap = (`mean_LLF_male'-`mean_LLF_female')-(`mean_NoLLF_male'-`mean_NoLLF_female')
	
	mat BF_effects_Third_Son = `change_in_gap'
	
	restore
	
}	

***********************************
* (3.0) Output table
***********************************
	putexcel B19 = "Implied Gender Gap in Breastfeeding, Late LLF Period"
		local iBFgap_1 = max(0, BF_effects_First[1,1])
		local iBFgap_1 = round(`iBFgap_1', .001)
		scalar cell = string(`iBFgap_1', "%05.3f")
		putexcel C19 = "`=cell'", hcenter

		local iBFgap_2s = max(0, BF_effects_Second_Son[1,1])
		local iBFgap_2s = round(`iBFgap_2s', .001)
		scalar cell = string(`iBFgap_2s', "%05.3f")
		putexcel D19 = "`=cell'", hcenter
		
		local iBFgap_2ns = max(0, BF_effects_Second_NoSon[1,1])
		local iBFgap_2ns = round(`iBFgap_2ns', .001)
		scalar cell = string(`iBFgap_2ns', "%05.3f")
		putexcel E19 = "`=cell'", hcenter
		
		local iBFgap_3s = max(0, BF_effects_Third_Son[1,1])
		local iBFgap_3s = round(`iBFgap_3s', .001)
		scalar cell = string(`iBFgap_3s', "%05.3f")
		putexcel F19 = "`=cell'", hcenter
		
		local iBFgap_3ns = max(0, BF_effects_Third_NoSon[1,1])
		local iBFgap_3ns = round(`iBFgap_3ns', .001)
		scalar cell = string(`iBFgap_3ns', "%05.3f")
		putexcel G19 = "`=cell'", hcenter

	putexcel B23 = "Breastfeeding Mortality Gap X Implied Gender Gap in Breastfeeding, Post LLF Period"
		local cell_1 = `mortGap_1'*`iBFgap_1'
		scalar cell_1 = string(`cell_1', "%05.3f")
		putexcel C23 = "`=cell_1'", hcenter
		
		local cell_2s = `mortGap_2s'*`iBFgap_2s'
		scalar cell = string(`cell_2s', "%05.3f")
		putexcel D23 = "`=cell'", hcenter
		
		local cell_2ns = `mortGap_2ns'*`iBFgap_2ns'
		scalar cell = string(`cell_2ns', "%05.3f")
		putexcel E23 = "`=cell'", hcenter
		
		local cell_3s = `mortGap_3s'*`iBFgap_3s'
		scalar cell = string(`cell_3s', "%05.3f")
		putexcel F23 = "`=cell'", hcenter
		
		local cell_3ns = `mortGap_3ns'*`iBFgap_3ns'
		scalar cell = string(`cell_3ns', "%05.3f")
		putexcel G23 = "`=cell'", hcenter
		
	putexcel B24 = "Excess Female Deaths Due to Implied Breastfeeing Gap"
		local excess_1 = round(`cell_1'*`births_1', 1)
		putexcel C24 = "`excess_1'", hcenter
		local excess_2s = round(`cell_2s'*`births_2s', 1)
		putexcel D24 = "`excess_2s'", hcenter
		local excess_2ns = round(`cell_2ns'*`births_2ns', 1)
		putexcel E24 = "`excess_2ns'", hcenter
		local excess_3s = round(`cell_3s'*`births_3s', 1)
		putexcel F24 = "`excess_3s'", hcenter
		local excess_3ns = round(`cell_3ns'*`births_3ns', 1)
		putexcel G24 = "`excess_3ns'", hcenter	

putexcel save

log close



